#Census data

#Define variables
agepovvars <- c("B01001_001", "B01001_020", "B01001_021", "B01001_022","B01001_023", "B01001_024","B01001_025", "B01001_044","B01001_045","B01001_046","B01001_047","B01001_048","B01001_049","B17001_015","B17001_016","B17001_029","B17001_030")

#grab data
ca <- get_acs(geography = "county",
              variables = agepovvars,
              state = "CA", 
              year = 2023)
## Getting data from the 2019-2023 5-year ACS
#widen stacked data
wide_ca <- ca %>%
  pivot_wider(names_from = variable, values_from = c(estimate, moe))

#If we wanted to eliminate the margins of error, we could do it all in a select statement:

pretty_ca <- wide_ca %>% dplyr::select(GEOID,NAME,pop=estimate_B01001_001, estimate_B01001_020, estimate_B01001_021, estimate_B01001_022, estimate_B01001_023, estimate_B01001_024, estimate_B01001_025, estimate_B01001_044, estimate_B01001_045, estimate_B01001_046, estimate_B01001_047, estimate_B01001_048,estimate_B01001_049,estimate_B17001_015,estimate_B17001_016,estimate_B17001_029,estimate_B17001_030)

#Calculate up_65, followed by percentages
ca_fin <- pretty_ca %>% rowwise() %>% mutate(up_65 = sum(estimate_B01001_020,estimate_B01001_021,estimate_B01001_022,estimate_B01001_023,estimate_B01001_024,estimate_B01001_025,estimate_B01001_044,estimate_B01001_045,estimate_B01001_046,estimate_B01001_047,estimate_B01001_048,estimate_B01001_049)) %>% 
mutate(pov_65=sum(estimate_B17001_015,estimate_B17001_016,estimate_B17001_029,estimate_B17001_030))

 
ca_fin <- ca_fin %>% mutate(pct_65up=100*(up_65/pop)) 
ca_fin <- ca_fin %>% mutate(pct_pov65=100*(pov_65/up_65))

#narrow the table to desired variables               
ca_fin <- ca_fin %>% dplyr::select(GEOID,NAME,pop,up_65,pov_65,pct_65up,pct_pov65)   
#grab map
cageo <- counties(state= 'CA', progress_bar = FALSE) 
## Retrieving data for the year 2021
#Join Data to map
camap <- merge(cageo, ca_fin, by.x = "GEOID", by.y = "GEOID")

#Do a qtm to check
qtm(camap,"pct_65up")

#Define a pop-up window:
capopup <- paste0("<b>COUNTY: ", camap$NAMELSAD, "</b><br />%65 and older: ", round(camap$pct_65up,2))

#Set color palette

capalette <- colorNumeric(palette = "Greens", domain=camap$pct_65up)

# Transform to leaflet projection if needed and map

camap <- st_transform(camap, crs = '+proj=longlat +datum=WGS84')

leaflet(camap) %>%
  addTiles() %>%
  addPolygons(stroke=FALSE, 
              smoothFactor = 0.2, 
              fillOpacity = .8, 
              popup=capopup, color= ~capalette(camap$pct_65up)) %>% 
              addLegend("bottomleft", pal = capalette, values = ~pct_65up,
    title = "<b>Percent 65 and older</b><br> ",
    opacity = 1
  )
#Define a pop-up window:
capopup2 <- paste0("<b>COUNTY: ", camap$NAMELSAD, "</b><br />%65 and older below poverty: ", round(camap$pct_pov65,2))

#Set color palette

capalette2 <- colorNumeric(palette = "Greens", domain=camap$pct_pov65)

# Transform to leaflet projection if needed and map

camap <- st_transform(camap, crs = '+proj=longlat +datum=WGS84')

leaflet(camap) %>%
  addTiles() %>%
  addPolygons(stroke=FALSE, 
              smoothFactor = 0.2, 
              fillOpacity = .8, 
              popup=capopup2, color= ~capalette2(camap$pct_pov65))%>% 
              addLegend("bottomleft", pal = capalette2, values = ~pct_pov65,
    title = "<b>Percent 65 and older below poverty</b><br> ",
    opacity = 1
  )